* VAR code for results in Leduc, S. and K. Sill, "Expectations and Economic Fluctuations: An Analysis Using Survey Data,"
* Review of Economic and Statistics 95(4), October 2013, 1352-1367.
* This code generates the responses to a positive shock to expected unemployment in the first column of Figure 10, which uses the Michigan expectation data.
* To reproduce the results in Figure 10, the positive expectation shock must be converted to a negative one.


********************************************************************************
*                                                                              *
* Note:  Computes 3 Confidence Intervals Simultaneously (90%, 80%, 68%)        *
*                                                                              *
********************************************************************************



********************************************************************************
* Note:  This program renormalizes the impulse responses to VAR shocks
*        by dividing all impulse responses by the own-contemporaneous response
*        (see the two "Patch Sections")
*
*        This program also computes historical decompositions & structural shock
*        estimates.
********************************************************************************

cal 1947 1 4
alloc 0 2011:1

****************************************
*                                      *
* Part 1. Read Data                    *
*                                      *
****************************************

open data c:\myfiles\papers\news\code\data_restat.rat
data(format=rats,org=obs) 1967:4 2010:4 umichm

****************************************************************
*                                                              *
* Part 2. Creating a new ordering for quarterly data           *
*                                                              *
****************************************************************

********************************
* a. First month of every quarter
********************************
data(format=rats,org=obs,select=1) 1947:1 2011:3 ruc rtb3 cpi rtb10
set ruc1      / = ruc(t)
set rtb31    / = rtb3(t)
set cpi1    / = cpi(t)
set rtb101     / = rtb10(t)

clear ruc rtb3 cpi rtb10

********************************
* a. Second month of every quarter
********************************
data(format=rats,org=obs,select=2) 1947:1 2011:3 ruc rtb3 cpi rtb10 umichm
set ruc2      / = ruc(t)
set rtb32    / = rtb3(t)
set cpi2    / = cpi(t)
set rtb102     / = rtb10(t)
set umichm2   / = umichm(t)

clear ruc rtb3 cpi rtb10 umichm

********************************
* a. Last month of every quarter
********************************
data(format=rats,org=obs,select=3) 1947:1 2011:3 ruc rtb3 cpi rtb10
set ruc3      / = ruc(t)
set rtb33     / = rtb3(t)
set cpi3     / = cpi(t)
set rtb103     / = rtb10(t)

clear ruc rtb3 cpi rtb10

********************************
* Quarterly data: Real-time RGDP
********************************
data(format=rats,org=obs) 1947:1 2011:3 grgdpunrevised


*************************************************************************************
*                                                                                   *
* Part 3. Compute VAR New Quarterly Macro Data                                      *
*                                                                                   *
*************************************************************************************


* Unemployment Rates
set urate   /  = (ruc2(t) + ruc3(t) + ruc1(t+1))/3.0

* Interest Rates
set r       /  = (rtb32(t) + rtb33(t) + rtb31(t+1))/3.0

* Inflation Rates
set inf         / = (log(cpi1(t+1)/cpi1(t)  ))*4.0*100.0

* 10-year T bill rate
set r10       /  = ((rtb102(t) + rtb103(t) + rtb101(t+1))/3.0)

* Real-time real GDP growth
set grgdprt / = grgdpunrevised(t)

* Oil Prices & Hamilton Dummy

*set oshock  /  = 0.0
*set oshock 1956:2 1956:2 = 10.1
*set oshock 1973:2 1973:2 =  7.8
*set oshock 1978:2 1978:2 =  8.9
*set oshock 1980:1 1980:1 =  7.2
*set oshock 1990:1 1990:1 =  8.8

*set oshockreplace  / = oshock

*set gshock  /  = 0.0
*set gshock 1980:1 1980:1 = 1


* Oil Shock Dummy
set oshock  /  = 0.0
set oshock 1973:4 1973:4 =  7.8
set oshock 1978:4 1978:4 =  8.9
set oshock 1980:4 1980:4 =  7.2
set oshock 1990:3 1990:3 =  8.8

set oshockreplace  / = oshock

* Fiscal Policy Shock Dummy
set gshock  /  = 0.0
set gshock 1980:1 1980:1 = 1
set gshock 2001:3 2001:3 = 1


************************************************************
*                                                          *
* Part 4.  Program Parameters For Lag Selection Routine    *
*                                                          *
************************************************************


*******************************************
* Set Generic Variables to Use in the VAR
*******************************************

set pdot      /  = inf;  				* actual inflation measure
set expec     /  = umichm2;				* expected unemployment rate
set u         /  = (ruc1(t)-ruc1(t-4)); 	* unemployment rate measure
set r         /  = r;  					* nominal rate measure


*********************************************************************************************
* Set Sample Periods & Other Parameters
*
* Note 1. We do not want lags of RHS variables to extend into the period before 1952:1
*         (in the first sample) or into the the period before 1979:2 (in the second sample).
*         Thus, given the samples listed below, we adjust the starting point of each sample
*         in the estimation commands, depending on the number of lags used in the VAR.
*********************************************************************************************

compute sbeg1 = 1978:1;
compute send1 = 2008:4;

declare rectangular[integer]  lagchoice(1,3)

        * Column 1  :   Corresponds to AIC (determined below)
        * Column 2  :   Corresponds to SIC (determined below)
        * Column 3  :   Any number you want

compute lagchoice(1,3) =  4 ; * enter an integer for number of lags to use, for use when you do
                              * not want to use the AIC or SIC (pre-Volcker period)

compute LS1   = 3;          ; * Lag Scheme(1 = AIC, 2 = SIC, 3 = other)

compute nvar  = 4           ; * number of variables in the VAR

                              * List order of variables as given below in the VARIABLES command

compute [vector[string]] varnames = ||'pdot', 'expec', 'u', 'r'||


******************************************************
*                                                    *
* Part 8.  Compute AIC & SIC Minimizing Lag Lengths  *
*                                                    *
******************************************************

set oshock1 / = oshock(t-1)
set oshock2 / = oshock(t-2)

*set gshock1 / = gshock(t-1)
*set gshock2 / = gshock(t-2)


*source(noecho) d:\myfiles\papers\news\code\aicsic.src


*@aicsic(more, nosuppress) sbeg1  send1  4  lagchoice(1,1) lagchoice(1,2)
*# pdot  expec u  r
*# constant oshock oshock1 oshock2
*# constant

*display;display;



                 ***************************************
                 *                                     *
                 *         Begin Bootstrap Routine     *
                 *                                     *
                 ***************************************



*********************************************************************************************
* This program produces the following:
*   IRPOINT_ = NVARS x NSHOCKS RECT[SERIES] of impulse response point estimates,
*             where the response of the ith variable to a shock to the jth equation is
*             IRPOINT(i,j)
*   LBAND_, UBAND_ =
*             NVARS x NSHOCKS RECT[SERIES] of lower and upper points of the confidence
*             interval, each dimensioned in the same way as IRPOINT
* Thus, the impulse response and confidence intervals for y(i) in response to a shock to
* equation j are given by IRPOINT_(i,j) and LBAND_(i,j), UBAND_(i,j).
*
* The (i,j) refer to the order in which the variables are listed in Part 4.a - 4.b (which is
* the same order that you must construct the VAR in Part 6), not the Cholesky ordering defined
* in Part 5.
*
*
*
* Note:  all bootstrap parameters and data types use "_" as the last character of the name
*
*
**********************************************************************************************



**********************************
*                                *
* Bootstap Parameters            *
*                                *
**********************************

compute sbeg1_     = sbeg1+lagchoice(1,LS1); * first period to use to estimate the VAR
                                             * (all lagged values MUST be defined at this date)

compute send1_     = send1;                  * last  period to use to estimate the VAR

compute subtitle_  = %datelabel(sbeg1_)+'-'+%datelabel(send1_)

compute ndraws_    = 1000;   * number of draws

compute nvars_     = nvar+1; * number of variables in the model = VAR eqns + identities
                             * (you must include at least one identity, thus nvars > nshocks)

compute nshocks_   = nvar;      * number of shocks in the VAR (i.e., nvars minus number of identities)
                             * (this is the number of equations in the VAR, not incl. identities)

compute nlags_     = lagchoice(1,LS1);   * number of lags in the VAR

compute nirsteps_  = 25;     * number of impulse response steps to compute

compute width1_     = 0.90;   * enter width for confidence interval 1 (0.95 means 95% intervals)
compute width2_     = 0.80;   * enter width for confidence interval 2 (0.95 means 95% intervals)
compute width3_     = 0.68;   * enter width for confidence interval 3 (0.95 means 95% intervals)


compute [vector] percentiles1_ $
                    = ||(1.0-width1_)/2.0,(1.0+width1_)/2.0||;  * percentiles for boostrapped error bands 1

compute [vector] percentiles2_ $
                    = ||(1.0-width2_)/2.0,(1.0+width2_)/2.0||;  * percentiles for boostrapped error bands 2

compute [vector] percentiles3_ $
                    = ||(1.0-width3_)/2.0,(1.0+width3_)/2.0||;  * percentiles for boostrapped error bands 3


**********************************
*                                *
* Dimension Objects              *
*                                *
**********************************


declare vector[series]        path_(nshocks_) actser_(nvars_)
declare rectangular           pathmatrix_(send1_-sbeg1_+1,nshocks_) irall_(ndraws_,nirsteps_*nvars_) $
                                                                    iroil_(ndraws_,nirsteps_*nvars_)

declare rectangular[series]   ssim_(ndraws_,nvars_) lband1_(nvars_,nshocks_) uband1_(nvars_,nshocks_)  $
                                                    lband2_(nvars_,nshocks_) uband2_(nvars_,nshocks_)  $
                                                    lband3_(nvars_,nshocks_) uband3_(nvars_,nshocks_)  $
                                                    irpointoil_(nvars_,1)                              $
                                                    lbandoil1_(nvars_,1)     ubandoil1_(nvars_,1)      $
                                                    lbandoil2_(nvars_,1)     ubandoil2_(nvars_,1)      $
                                                    lbandoil3_(nvars_,1)     ubandoil3_(nvars_,1)



declare vector[integer]       sbegs_(nvars_)




*****************************************************************
*                                                               *
* Place All Relevant Data Into VEC[SER] = ACTSER_               *
*                                                               *
*   Note 1. The ordering does not affect the Cholesky Decomp    *
*   Note 2. The ordering only affects the position of the VAR   *
*           equations in the model                              *
*   Note 3. You must have nvars entries, VAR equation variables *
*           first, then variables defined by identity           *
*                                                               *
*****************************************************************

* 4.A.  VAR Equation Variables (these are the dependent variables of the VAR)

set actser_(1)    /  = pdot
set actser_(2)    /  = expec
set actser_(3)    /  = u
set actser_(4)    /  = r



do i = 1, nshocks_
   labels actser_(i)
   # varnames(i)
end do i


* 4.B.  Variables Defined By Identity (these are the dependent variables in the identities)

set realrate      /  = r(t) - pdot(t)

set actser_(5)    /  = realrate(t)

* 4.C. Create Vector Of Starting Dates For Each ACTSER_(i)

do i = 1,nvars_
   inquire(series=actser_(i)) sbegs_(i) scrap
end do i


***************************************
*                                     *
* Define Cholesky Ordering            *
*                                     *
***************************************

compute [vector[integer]] corder_ = ||2,1,4,3||
compute choleskystring_ = ' / Order = U, Pdot, R, expec'



**********************************************************
*                                                        *
* Define The VAR                                         *
*                                                        *
*    Note 1. List variables in the order given in 4.A    *
*                                                        *
**********************************************************


system(model=VAR1eqns_)
     variables  actser_(1) actser_(2) actser_(3) actser_(4)
     lags 1 to nlags_
     *det constant oshock{0 to nlags_}
     det constant oshock{0 to nlags_} gshock{0 to nlags_}
      *det constant
end(system)



***********************************************
*                                             *
* Define All Identities                       *
*                                             *
*                                             *
* Note: you must hace at least one identity   *
***********************************************



* Real Rate
equation(identity) realrate_id   actser_(5)
# actser_(4)  actser_(1)
associate          realrate_id
# 1.0        -1.0



*********************************************************************
*                                                                   *
* Estimate The VAR Model                                            *
*         Capture Residuals                                         *
*         Capture Coefficients                                      *
*         Compute Cholesky Factor                                   *
*         Group Estimated VAR Equations & Identities Into a Model   *
*                                                                   *
*********************************************************************

estimate(residuals=resid_,coeffs=beta_)  sbeg1_  send1_
*linreg expec sbeg1_ send1_ resids
*# constant expec{1} expec{2} pdot{1} pdot{2} u{1} u{2} r{1} r{2} oshock{0 to 2} gshock{0 to 2}
*statistics resids

*linreg expec sbeg1_ send1_
*# constant expec{1 to 2} pdot{1 to 2} u{1 to 2} r{1 to 2} oshock{0 to 2} gshock{0 to 2}
*exclude(title="Granger Causality Test")
*# pdot{1 to 2} u{1 to 2} r{1 to 2} oshock{0 to 2} gshock{0 to 2}


compute cfactor_ = %psdfactor(%sigma,corder_)


group ID1model_  realrate_id

compute [model] model1_     = VAR1eqns_ + ID1model_


*******************************************************************************
*                                                                             *
* Compute & Capture Impulse Response Point Estimates, Variance Decomps,       *
* Historical Decompositions & Structural Shocks                               *
*                                                                             *
* Note 1.  The impulse response point estimates are in the NVARS_ x NSHOCKS_  *
*          RECT[SERIES] called IRPOINT_(i,j), i = 1,...,NVARS_ and            *
*          j = 1,...,NSHOCKS_.  The response of variable i to a shock         *
*          to equation j is IRPOINT_(i,j), where the number is that given in  *
*          Part 4A-4B, not by the Cholesky chosen in Part 5                   *
*                                                                             *
*******************************************************************************

display 'Cholesky Factor = ' cfactor_

****************************************
* Compute Structural Errors
*    = S_ERROR, a nshocks_ x 1 VEC[SER]
****************************************

declare vector[integer]  residlist
declare vector[series]   s_error(nshocks_)

enter(varying) residlist
# resid_(1)
do i = 2, nshocks_
   enter(varying) residlist
   # residlist resid_(i)
end do i

make rf_errormat  sbeg1_  send1_
# residlist

compute s_errormat = (rf_errormat)*(tr(inv(cfactor_)))

do i = 1, nshocks_
   set s_error(i) sbeg1_ send1_ = s_errormat(t-(sbeg1_-1),i)
end do i



*******************************************************************
* Compute: Var Decomp, Hist Decomp, Impulse Responses
* Note: HISTORY command will not work on a model with identities
*******************************************************************

errors(model =model1_,decomp=cfactor_)                      *  nirsteps_  *
history(model=VAR1eqns_,decomp=cfactor_,results=hist_,noadd)  *  (send1_-sbeg1_+1) sbeg1_
impulse(model=model1_,decomp=cfactor_,results=irpoint_)     *  nirsteps_  *

print(dates) / hist_

******************************************************************************************
*         PATCH:  divide all impulse responses by contemporaneous own response
******************************************************************************************
do j = 1, nshocks_
   compute rescale = irpoint_(j,j)(1)
   do i = 1, nvars_
      set irpoint_(i,j) 1 nirsteps_ = irpoint_(i,j)/rescale
   end do i
end do j
******************************************************************************************



*******************************************************************************
*                                                                             *
* Compute & Capture Impulse Response Point Estimates to Oil Dummy Shock       *
*                                                                             *
* Note:  we compute this as the difference between 2 forecasts, one of which  *
*        has a (temporary) unit increase in Hamilton's oil dummy measure in   *
*        the first period of the forecast.                                    *
*                                                                             *
*******************************************************************************

****************************
* Compute Baseline Forecast
****************************

forecast(model=model1_,results = noshockfor_) * nirsteps_  sbeg1_

***********************************************
* Compute Forecast With Higher Dummy Variable
***********************************************

set oshock sbeg1_ sbeg1_ = oshock + 1.0

forecast(model=model1_,results = shockfor_)   * nirsteps_  sbeg1_

set oshock sbeg1_ sbeg1_ = oshockreplace


*************************************
* Compute Implied Impulse Responses (as the difference between the forecasts)
*************************************

do i = 1, nvars_

   set irpointoil_(i,1) / = shockfor_(i) - noshockfor_(i)

end do i

clear shockfor_  noshockfor_



*******************************************************************************
*                                                                             *
* Compute Simulated VAR Data                                                  *
*                                                                             *
* Note 1. The simulated values are in the NDRAWS x NVARS RECT[SERIES]         *
*         called SSIM_(draw,i),where the row (draw) indicates the replication *
*         number and the column (i) indicates the variable simulated          *
*                                                                             *
*******************************************************************************

* Conduct Simulation

smpl sbeg1_ send1_
do draw_ = 1, ndraws_

   boot entries_ / sbeg1_ send1_

   do i = 1,nshocks_
      set path_(i) / = resid_(i)(entries_(t))
   end do i

   do i = 1, nshocks_
      do j = 1, send1_-sbeg1_+1
         compute pathmatrix_(j,i) = path_(i)(sbeg1_-1+j)
      end do j
   end do i

   forecast(model=model1_,matrix=pathmatrix_, results=results_)

   do i = 1,nvars_
      set ssim_(draw_,i) / = results_(i)
   end do i

end do draw_


smpl
* Reset SMPL Above This Line


*********************************************************************************************
*                                                                                           *
* Estimate VAR  on Simulated Data                                                           *
*          Compute Impulse Responses                                                        *
*          Calculate CIs                                                                    *
*                                                                                           *
* Note 1.  The lower and upper confidence intervals are given in the NVARS x NSHOCKS        *
*          RECT[SERIES] called LBAND_(var,shock) and UBAND_(var,shock), where the row (var) *
*          indicates the response variable and the column (shock) indicates the equation    *
*          that is shocked. Note that the variable number and equation-shock number are     *
*          given in Part 4.A-4.B, not the Cholesky ordering in Part 5.                      *
*                                                                                           *
*********************************************************************************************

* Splice Initial Conditions (given by history) On Top Of Simulated Values

do draw_ = 1, ndraws_
   do var_ = 1, nvars_

      set ssim_(draw_,var_) sbegs_(var_) sbeg1_-1  = actser_(var_)

   end do var_
end do draw_


* Impulse Responses

do shock_ = 1,nshocks_

do draw_=1,ndraws_

   * Estimation

   system(model=VAR2eqns_)
     variables ssim_(draw_,1) ssim_(draw_,2) ssim_(draw_,3) ssim_(draw_,4)
     lags 1 to nlags_
     *det constant oshock{0 to nlags_}
     det constant oshock{0 to nlags_} gshock{0 to nlags_}
     *det constant
   end(system)

   estimate(noprint,noftests,coeffs=simbeta_) sbeg1_ send1_



   * Compute Coefficient Averages (over draws)[Not Used In This Program]

   if draw_.eq.1
     {
       compute musimbeta_ = (1.0/ndraws_)*simbeta_
     }
   else
     {
       compute musimbeta_ = musimbeta_ + (1.0/ndraws_)*simbeta_
     }

   * Define Identities

   * Real Rate
   equation(identity) realrate_id   ssim_(draw_,5)
   # ssim_(draw_,4)  ssim_(draw_,1)
   associate          realrate_id
   # 1.0           -1.0


   group ID2model_  realrate_id

   compute [model] model2_     = VAR2eqns_ + ID2model_


  * Compute Impulse Responses

   * Note 1: the equation shocked by the loop index, called "shock", is given by the order of
   *         equations in the VAR, not the order given by the cholesky factorization
   * Note 2: IR is an Nvars x 1 RECTANG[SER], with the ordering given by the ordering of the
   *         equations in the VAR, not the order of the cholesky factorization, and the
   *         ordering of the identities


   compute cfactor_ = %psdfactor(%sigma,corder_)
   impulse(model=model2_,noprint,decomp=cfactor_,results=ir_) * nirsteps_ shock_




   *****************************************************************************************
   *        PATCH:  divide all impulse responses by contemporaneous own response
   *****************************************************************************************
   compute rescale = ir_(shock_,1)(1)
   do i = 1, nvars_
     set ir_(i,1) 1 nirsteps_ = ir_(i,1)/rescale
   end do i
   *****************************************************************************************




   * Capture Impulse Responses For This Draw


   do var_ = 1,nvars_
      do step_ = 1,nirsteps_
          compute irall_(draw_,step_+(var_-1)*nirsteps_) = ir_(var_,1)(step_)
      end do step_
   end do var_

  * Compute & Capture Impulse Responses For Oil Dummy

   if shock_.eq.1
   {

      * Compute Baseline Forecast

      forecast(model=model2_,results = noshockfor_) * nirsteps_  sbeg1_

      * Compute Forecast With Higher Dummy Variable

      set oshock sbeg1_ sbeg1_ = oshock + 1.0

      forecast(model=model2_,results = shockfor_)   * nirsteps_  sbeg1_

      set oshock sbeg1_ sbeg1_ = oshockreplace

      * Compute Impulse Response As Difference Between The Forecasts & Capture

      do var_ = 1,nvars_
         do step_ = 1,nirsteps_
            compute iroil_(draw_,step_+(var_-1)*nirsteps_) = shockfor_(var_)(sbeg1_-1+step_)   $
                                                           - noshockfor_(var_)(sbeg1_-1+step_)
         end do step_
      end do var_

      clear noshockfor_ shockfor_

   }

end do draw_


   * Compute Confidence Intervals For This VAR Shock & The Oil Dummy Shock (by variable and step)

   do var_ = 1,nvars_
      do step_ = 1,nirsteps_

         * VAR Shock Responses
         compute col_ = %xcol(irall_,step_+(var_-1)*nirsteps_)
         compute lowerupper1_ = %fractiles(col_,percentiles1_)
         compute lowerupper2_ = %fractiles(col_,percentiles2_)
         compute lowerupper3_ = %fractiles(col_,percentiles3_)


         set lband1_(var_,shock_) step_ step_ = lowerupper1_(1)
         set lband2_(var_,shock_) step_ step_ = lowerupper2_(1)
         set lband3_(var_,shock_) step_ step_ = lowerupper3_(1)


         set uband1_(var_,shock_) step_ step_ = lowerupper1_(2)
         set uband2_(var_,shock_) step_ step_ = lowerupper2_(2)
         set uband3_(var_,shock_) step_ step_ = lowerupper3_(2)


         * Oil Dummy Shock Responses
         if shock_.eq.1
         {
            compute coloil_ = %xcol(iroil_,step_+(var_-1)*nirsteps_)

            compute lowerupperoil1_ = %fractiles(coloil_,percentiles1_)
            compute lowerupperoil2_ = %fractiles(coloil_,percentiles2_)
            compute lowerupperoil3_ = %fractiles(coloil_,percentiles3_)

            set lbandoil1_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil1_(1)
            set lbandoil2_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil2_(1)
            set lbandoil3_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil3_(1)

            set ubandoil1_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil1_(2)
            set ubandoil2_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil2_(2)
            set ubandoil3_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil3_(2)

         }

      end do step_
   end do var_

end do shock_



*********************************************
*                                           *
* Graphs of Impulse Responses:  VAR Shocks  *
*                                           *
*********************************************


*********************
* Expectations Shock
*********************

spgraph(vfields=3,hfields=2,header='Expectations Shock')

   graph(header='Inflation Response',nodates,number=0) 7
   # lband1_(1,2)     1 nirsteps_   1
   # irpoint_(1,2)    1 nirsteps_   4
   # uband1_(1,2)     1 nirsteps_   1
   # lband2_(1,2)     1 nirsteps_   3
   # uband2_(1,2)     1 nirsteps_   3
   # lband3_(1,2)     1 nirsteps_   6
   # uband3_(1,2)     1 nirsteps_   6

   graph(header='Expectations Response',nodates,number=0) 7
   # lband1_(2,2)     1 nirsteps_   1
   # irpoint_(2,2)    1 nirsteps_   4
   # uband1_(2,2)     1 nirsteps_   1
   # lband2_(2,2)     1 nirsteps_   3
   # uband2_(2,2)     1 nirsteps_   3
   # lband3_(2,2)     1 nirsteps_   6
   # uband3_(2,2)     1 nirsteps_   6

   graph(header='Unemployment Response',nodates,number=0) 7
   # lband1_(3,2)     1 nirsteps_   1
   # irpoint_(3,2)    1 nirsteps_   4
   # uband1_(3,2)     1 nirsteps_   1
   # lband2_(3,2)     1 nirsteps_   3
   # uband2_(3,2)     1 nirsteps_   3
   # lband3_(3,2)     1 nirsteps_   6
   # uband3_(3,2)     1 nirsteps_   6

   graph(header='3-month Interest Rate Response',nodates,number=0) 7
   # lband1_(4,2)     1 nirsteps_   1
   # irpoint_(4,2)    1 nirsteps_   4
   # uband1_(4,2)     1 nirsteps_   1
   # lband2_(4,2)     1 nirsteps_   3
   # uband2_(4,2)     1 nirsteps_   3
   # lband3_(4,2)     1 nirsteps_   6
   # uband3_(4,2)     1 nirsteps_   6

spgraph(done)


*******************************************
*                                         *
* Graphs of Historical Decompositions     *
*                                         *
*******************************************

* Note:  The decomposition for variable i is in column i.  The first row is the base projection.
*        Row 2: component due to first variable in the VAR (pdot)
*        Row 3: component due to second variable in the VAR (expec), etc.



*************************************
* Inflation Historical Decomposition
*************************************

*set error_pdot = pdot - hist_(1,1)

*spgraph(vfields=3,hfields=2,header='Cumulative Effects of Shocks on Inflation', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Actual Inflation (solid) & Base Projection (dash)',patterns) 2
*    # pdot         sbeg1_  send1_   1
*    # hist_(1,1)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Inflation Shocks (dash)',patterns)    2
*    # error_pdot   sbeg1_  send1_   1
*    # hist_(2,1)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Expectations Shocks (dash)',patterns) 2
*    # error_pdot   sbeg1_  send1_   1
*    # hist_(3,1)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Unemployment Shocks (dash)',patterns) 2
*    # error_pdot   sbeg1_  send1_   1
*    # hist_(4,1)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Monetary Policy Shocks (dash)',patterns) 2
*    # error_pdot   sbeg1_  send1_   1
*    # hist_(5,1)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Commodity Price Shocks (dash)',patterns) 2
*    # error_pdot   sbeg1_  send1_   1
*    # hist_(6,1)   sbeg1_  send1_   2

*spgraph(done)


****************************************
* Expectations Historical Decomposition
****************************************

*set error_expec = expec - hist_(1,2)

*spgraph(vfields=3,hfields=2,header='Cumulative Effects of Shocks on Expectations', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Actual Expectations (solid) & Base Projection (dash)',patterns) 2
*    # expec        sbeg1_  send1_   1
*    # hist_(1,2)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Inflation Shocks (dash)',patterns)    2
*    # error_expec  sbeg1_  send1_   1
*    # hist_(2,2)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Expectations Shocks (dash)',patterns) 2
*    # error_expec  sbeg1_  send1_   1
*    # hist_(3,2)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Unemployment Shocks (dash)',patterns) 2
*    # error_expec  sbeg1_  send1_   1
*    # hist_(4,2)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Monetary Policy Shocks (dash)',patterns) 2
*    # error_expec  sbeg1_  send1_   1
*    # hist_(5,2)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Commodity Price Shocks (dash)',patterns) 2
*    # error_expec  sbeg1_  send1_   1
*    # hist_(6,2)   sbeg1_  send1_   2

*spgraph(done)


****************************************
* Unemployment Historical Decomposition
****************************************

*set error_u = u - hist_(1,3)

*spgraph(vfields=3,hfields=2,header='Cumulative Effects of Shocks on Unemployment', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Actual Unemployment(solid) & Base Projection (dash)',patterns) 2
*    # u            sbeg1_  send1_   1
*    # hist_(1,3)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Inflation Shocks (dash)',patterns)    2
*    # error_u      sbeg1_  send1_   1
*    # hist_(2,3)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Expectations Shocks (dash)',patterns) 2
*    # error_u      sbeg1_  send1_   1
*    # hist_(3,3)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Unemployment Shocks (dash)',patterns) 2
*    # error_u      sbeg1_  send1_   1
*    # hist_(4,3)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Monetary Policy Shocks (dash)',patterns) 2
*    # error_u      sbeg1_  send1_   1
*    # hist_(5,3)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Commodity Price Shocks (dash)',patterns) 2
*    # error_u      sbeg1_  send1_   1
*    # hist_(6,3)   sbeg1_  send1_   2

*spgraph(done)


*****************************************
* Nominal Rate Historical Decomposition
*****************************************

*set error_r = r - hist_(1,4)

*spgraph(vfields=3,hfields=2,header='Cumulative Effects of Shocks on the Nominal Rate', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Actual Rate (solid) & Base Projection (dash)',patterns) 2
*    # r            sbeg1_  send1_   1
*    # hist_(1,4)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Inflation Shocks (dash)',patterns)    2
*    # error_r      sbeg1_  send1_   1
*    # hist_(2,4)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Expectations Shocks (dash)',patterns) 2
*    # error_r      sbeg1_  send1_   1
*    # hist_(3,4)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Unemployment Shocks (dash)',patterns) 2
*    # error_r      sbeg1_  send1_   1
*    # hist_(4,4)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Monetary Policy Shocks (dash)',patterns) 2
*    # error_r      sbeg1_  send1_   1
*    # hist_(5,4)   sbeg1_  send1_   2

*   graph(header='Deviation From Base (line) & Effects of Commodity Price Shocks (dash)',patterns) 2
*    # error_r      sbeg1_  send1_   1
*    # hist_(6,4)   sbeg1_  send1_   2

*spgraph(done)



*******************************************
* Commodity Price Historical Decomposition
*******************************************

*set error_pcom = pcom - hist_(1,5)

*spgraph(vfields=3,hfields=2,header='Cumulative Effects of Shocks on Commodity Prices', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Actual Rate (solid) & Base Projection (dash)',patterns) 2
*    # pcom         sbeg1_  send1_   1
*    # hist_(1,5)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Inflation Shocks (dash)',patterns)    2
*    # error_pcom   sbeg1_  send1_   1
*    # hist_(2,5)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Expectations Shocks (dash)',patterns) 2
*    # error_pcom   sbeg1_  send1_   1
*    # hist_(3,5)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Unemployment Shocks (dash)',patterns) 2
*    # error_pcom   sbeg1_  send1_   1
*    # hist_(4,5)   sbeg1_  send1_   2

*    graph(header='Deviation From Base (line) & Effects of Monetary Policy Shocks (dash)',patterns) 2
*    # error_pcom   sbeg1_  send1_   1
*    # hist_(5,5)   sbeg1_  send1_   2

*   graph(header='Deviation From Base (line) & Effects of Commodity Price Shocks (dash)',patterns) 2
*    # error_pcom   sbeg1_  send1_   1
*    # hist_(6,5)   sbeg1_  send1_   2

*spgraph(done)




***********************************
*                                 *
* Graphs of Structural Shocks     *
*                                 *
***********************************


*spgraph(vfields=3,hfields=2,header='Estimates of Structural Shocks', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Inflation Shocks')       1
*    # s_error(1)   sbeg1_  send1_ 1

*    graph(header='Expectations Shocks')    1
*    # s_error(2)   sbeg1_  send1_ 1

*    graph(header='Unemployment Shocks')    1
*    # s_error(3)   sbeg1_  send1_ 1

*    graph(header='Monetary Policy Shocks') 1
*    # s_error(4)   sbeg1_  send1_ 1

*    graph(header='Commodity Price Shocks') 1
*    # s_error(5)   sbeg1_  send1_ 1

*spgraph(done)





**************************************************************
*                                                            *
* Output Impulse Responses & Hist Decomps to a Worksheet     *
*                                                            *
**************************************************************

*******************************
* VAR Shock Impulse Responses
*******************************

*set steps 1 nirsteps_ = t-1

* Confidence Intervals 1

*open copy c:\myfiles\papers\news\code\output\March08\ordered_first\new_timing\6months_ahead\dec09\revisions\VarIR_MICH4WDUMMY_BaseO_90.xls
*copy(nodates,format=xls,org=obs) 1  nirsteps_  $
*                                 steps  irpoint_  lband1_   uband1_

* Confidence Intevals 2

*open copy d:\myfiles\papers\sunspots\Tom's code\output\VarIR_PreV_Prg7B1_BaseO_80.xls
*copy(nodates,format=xls,org=obs) 1  nirsteps_  $
*                                 steps  irpoint_  lband2_   uband2_

* Confidence Intervals 3

*open copy c:\myfiles\papers\news\code\output\March08\ordered_first\new_timing\6months_ahead\dec09\revisions\VarIR_MICH4WDUMMY_BaseO_68.xls
*copy(nodates,format=xls,org=obs) 1  nirsteps_  $
*                                 steps  irpoint_  lband3_   uband3_



***********************
* Historical Decomps
***********************

*open copy c:\myfiles\papers\sunspots\Tom's code\output\Hist_PreV_Prg7B1_BaseO.xls
*copy(dates,format=xls,org=obs) /                           $
*                                  pdot  expec  u  r  pcom  $
*                                  hist_



























































































































































